7-2 WHO 8-point Outcome Scale

Analyses of the WHO 8-point outcome scale outcome.

Authors

Affiliations

James Totterdell

University of Sydney

Rob Mahar

University of Melbourne

Load packages

library(ASCOTr)
library(tidyverse)
library(lubridate)
library(kableExtra)
library(patchwork)
library(cmdstanr)
library(posterior)
library(bayesplot)
library(ggdist)
library(lme4)
library(broom)
library(broom.mixed)
library(bayestestR)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))

bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))

color_scheme_set("red")
options(digits = 4)

Prepare datasets

all_data <- readRDS(file.path(ASCOT_DATA, "all_data.rds"))

# Create censored D28 outcome for the two participants
all_data <- all_data |> mutate(
  D28_who_low = case_when(
    StudyPatientID %in% c("ALF00006", "ALF00012") ~ 1L,
    TRUE ~ D28_who
  ),
  D28_who_high = case_when(
    StudyPatientID %in% c("ALF00006", "ALF00012") ~ 2L,
    TRUE ~ D28_who
  )
)

# FAS-ITT
fas_itt_dat <- ASCOTr:::make_fas_itt_set(all_data)
fas_itt_nona_dat <- fas_itt_dat |>
  filter(!is.na(D28_who))

# ACS-ITT
acs_itt_dat <- ASCOTr:::make_acs_itt_set(all_data)
acs_itt_nona_dat <- acs_itt_dat |>
  filter(!is.na(D28_who))

# AVS-ITT
avs_itt_dat <- ASCOTr:::make_avs_itt_set(all_data)
avs_itt_nona_dat <- avs_itt_dat |>
  filter(!is.na(D28_who))

Load models

ordmod0 <- compile_cmdstanr_mod(
  file.path("ordinal", "logistic_cumulative"), dir = "stan")
ordmodcens0 <- compile_cmdstanr_mod(
  file.path("ordinal", "logit_cumulative_cens"), dir = "stan")
ordmod_site <- compile_cmdstanr_mod(
  file.path("ordinal", "logistic_cumulative_site"), dir = "stan")
ordmod <- compile_cmdstanr_mod(
  file.path("ordinal", "logistic_cumulative_epoch_site"), dir = "stan")
ordmodcens <- compile_cmdstanr_mod(
  file.path("ordinal", "logit_cumulative_epoch_site_cens"), dir = "stan")
logistic <- compile_cmdstanr_mod(
  file.path("binary", "logistic_site_epoch"), dir = "stan")

Helper functions

 make_primary_model_data2 <- function(dat,
                                    outcome = "PO",
                                    vars = c("inelgc3", "agegte60", "ctry"),
                                    beta_sd_int = 2.5,
                                    beta_sd_var = c(10, 2.5, 1, 1),
                                    beta_sd_trt = 1,
                                    ctr = contr.equalprior,
                                    includeA = TRUE,
                                    includeC = TRUE,
                                    intercept = TRUE,
                                    ...) {
  X <- make_X_design(dat, vars = vars, ctr = ctr, includeA = includeA, includeC = includeC, intercept = intercept)
  nXtrt <- sum(grepl("rand", colnames(X)))
  epoch <- dat[["epoch"]]
  M_epoch <- max(epoch)
  region <- dat[["ctry_num"]]
  M_region <- max(region)
  site <- dat[["site_num"]]
  M_site <- max(site)
  region_by_site <- region_by_site <- dat |>
    dplyr::count(ctry_num, site_num) |>
    pull(ctry_num)
  y_raw <- as.matrix(dat[, outcome])
  y_mod <- y_raw
  N <- dim(X)[1]
  K <- dim(X)[2]
  beta_sd <- c(rep(beta_sd_trt, nXtrt), beta_sd_var)
  if (intercept) {
    beta_sd <- c(beta_sd_int, beta_sd)
  }
  out <- list(
    N = N, K = K, X = X, y = y_mod, y_raw = y_raw,
    J = max(y_mod), p_par = rep(2 / max(y_mod), max(y_mod)),
    M_region = M_region, region = region,
    M_site = M_site, site = site,
    M_epoch = M_epoch, epoch = epoch,
    region_by_site = region_by_site,
    beta_sd = beta_sd
  )
  if (includeA) {
    out <- c(out, list(ctrA = attr(X, "contrasts")[["randA"]]))
  }
  if (includeC) {
    out <- c(out, list(ctrC = attr(X, "contrasts")[["randC"]]))
  }
  return(out)
 }

fit_primary_model2 <- function(dat = NULL,
                              model = NULL,
                              outcome = "PO",
                              vars = c("inelgc3", "agegte60", "sexF", "supp_oxy2", "ctry"),
                              beta_sd_int = 2.5,
                              beta_sd_var = c(10, 2.5, 2.5, 2.5, 1, 1),
                              beta_sd_trt = 1,
                              ctr = contr.equalprior,
                              includeA = TRUE,
                              includeC = TRUE,
                              intercept = TRUE,
                              seed = 32915,
                              adapt_delta = 0.99,
                              iter_sampling = 2500,
                              chains = 8,
                              ...) {
  mdat <- make_primary_model_data2(
    dat,
    outcome = outcome,
    vars = vars,
    beta_sd_int = beta_sd_int,
    beta_sd_var = beta_sd_var,
    beta_sd_trt = beta_sd_trt,
    ctr = ctr,
    includeA = includeA,
    includeC = includeC,
    intercept = intercept,
    ...
  )
  snk <- capture.output(
    mfit <- model[["sample"]](
      data = mdat,
      seed = seed,
      adapt_delta = adapt_delta,
      refresh = 0,
      iter_warmup = 1000,
      iter_sampling = iter_sampling,
      chains = chains,
      parallel_chains = min(chains, parallel::detectCores())
    )
  )
  mpars <- mfit$metadata()$model_params
  keep <- mpars[!grepl("(_raw|epsilon_)", mpars)]
  mdrws <- as_draws_rvars(mfit$draws(keep))
  names(mdrws$beta) <- colnames(mdat$X)
  if (any(grepl("site", keep))) {
    site_map <- dat %>% dplyr::count(site_num, site)
    names(mdrws$gamma_site) <- site_map$site
  }
  if (any(grepl("epoch", keep))) {
    epoch_map <- dat %>% dplyr::count(epoch, epoch_lab)
    names(mdrws$gamma_epoch) <- epoch_map$epoch_lab
  }

  if (includeA) {
    Ca <- mdat$ctrA
    mdrws$Acon <- as.vector(Ca %**% mdrws$beta[grepl("randA[0-9]", names(mdrws$beta))])
    names(mdrws$Acon) <- rownames(Ca)
    mdrws$Atrt <- mdrws$Acon[-1] - mdrws$Acon[1]
    mdrws$AOR <- exp(mdrws$Atrt)
  }
  if (includeC) {
    Cc <- mdat$ctrC
    mdrws$Ccon <- as.vector(Cc %**% mdrws$beta[grepl("randC[0-9]", names(mdrws$beta))])
    names(mdrws$Ccon) <- rownames(Cc)
    mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
    mdrws$COR <- exp(mdrws$Ctrt)
  }
  mdrws$OR <- exp(mdrws$beta[!grepl("(Intercept|rand)", names(mdrws$beta))])
  return(list(dat = mdat, fit = mfit, drws = mdrws))
}

make_summary_table_anticoagulation <- function(dat, format = "html") {
  tdat <- dat %>%
  group_by(
    CAssignment = factor(CAssignment, labels = intervention_labels2()$CAssignment)) %>%
  summarise(
    Patients = n(),
    Known = sum(!is.na(D28_who)),
    Deaths = sprintf("%i (%.0f%%)", 
                     sum(D28_who == 8, na.rm = TRUE), 
                     100 * mean(D28_who == 8, na.rm = TRUE)),
    `Hospitalised` = sprintf("%i (%.0f%%)", 
                     sum(D28_who > 2 & D28_who < 8, na.rm = TRUE), 
                     100 * mean(D28_who > 2 & D28_who < 8, na.rm = TRUE)),
    `WHO, Median (Q1, Q3)` = sprintf(
      "%.0f (%.0f, %.0f)", 
      median(D28_who, na.rm = T), 
      quantile(D28_who, 0.25, na.rm = TRUE), 
      quantile(D28_who, 0.75, na.rm = TRUE))
  ) %>%
  bind_rows(
    dat %>%
  group_by(CAssignment = "Overall") %>%
  summarise(
    Patients = n(),
    Known = sum(!is.na(D28_who)),
    Deaths = sprintf("%i (%.0f%%)", 
                     sum(D28_who == 8, na.rm = TRUE), 
                     100 * mean(D28_who == 8, na.rm = TRUE)),
    `Hospitalised` = sprintf("%i (%.0f%%)", 
                     sum(D28_who > 2 & D28_who < 8, na.rm = TRUE), 
                     100 * mean(D28_who > 2 & D28_who < 8, na.rm = TRUE)),
    `WHO, Median (Q1, Q3)` = sprintf(
      "%.0f (%.0f, %.0f)", 
      median(D28_who, na.rm = T), 
      quantile(D28_who, 0.25, na.rm = TRUE), 
      quantile(D28_who, 0.75, na.rm = TRUE))
  )
  ) %>%
  rename(`Anticoagulation\nintervention` = CAssignment)
  kable(
    tdat,
    format = format,
    align = "lrrrrr",
    booktabs = TRUE,
    linesep = ""
  ) %>%
    kable_styling(
      font_size = 9,
      latex_options = "HOLD_position"
    ) %>%
    row_spec(nrow(tdat), bold = T)
}

make_summary_table_antiviral <- function(dat, format = "html") {
  tdat <- dat %>%
  group_by(
    AAssignment = factor(AAssignment, labels = intervention_labels2()$AAssignment)) %>%
  summarise(
    Patients = n(),
    Known = sum(!is.na(D28_who)),
    Deaths = sprintf("%i (%.0f%%)", 
                     sum(D28_who == 8, na.rm = TRUE), 
                     100 * mean(D28_who == 8, na.rm = TRUE)),
    `Hospitalised` = sprintf("%i (%.0f%%)", 
                     sum(D28_who > 2 & D28_who < 8, na.rm = TRUE), 
                     100 * mean(D28_who > 2 & D28_who < 8, na.rm = TRUE)),
    `WHO, Median (Q1, Q3)` = sprintf(
      "%.0f (%.0f, %.0f)", 
      median(D28_who, na.rm = T), 
      quantile(D28_who, 0.25, na.rm = TRUE), 
      quantile(D28_who, 0.75, na.rm = TRUE))
  ) %>%
  bind_rows(
    dat %>%
  group_by(AAssignment = "Overall") %>%
  summarise(
    Patients = n(),
    Known = sum(!is.na(D28_who)),
    Deaths = sprintf("%i (%.0f%%)", 
                     sum(D28_who == 8, na.rm = TRUE), 
                     100 * mean(D28_who == 8, na.rm = TRUE)),
    `Hospitalised` = sprintf("%i (%.0f%%)", 
                     sum(D28_who > 2 & D28_who < 8, na.rm = TRUE), 
                     100 * mean(D28_who > 2 & D28_who < 8, na.rm = TRUE)),
    `WHO, Median (Q1, Q3)` = sprintf(
      "%.0f (%.0f, %.0f)", 
      median(D28_who, na.rm = T), 
      quantile(D28_who, 0.25, na.rm = TRUE), 
      quantile(D28_who, 0.75, na.rm = TRUE))
  )
  ) %>%
  rename(`Antiviral\nintervention` = AAssignment)
  kable(
    tdat,
    format = format,
    align = "lrrrrr",
    booktabs = TRUE,
    linesep = ""
  ) %>%
    kable_styling(
      font_size = 9,
      latex_options = "HOLD_position"
    ) %>%
    row_spec(nrow(tdat), bold = T)
}

Descriptive

Anticoagulation

Summary of WHO outcome by arm

save_tex_table(
  make_summary_table_anticoagulation(fas_itt_dat, "latex"), 
  file.path("outcomes", "secondary", "7-2-anticoagulation-summary"))
make_summary_table_anticoagulation(fas_itt_dat)

Anticoagulation intervention	Patients	Known	Deaths	Hospitalised	WHO, Median (Q1, Q3)
Not randomised to anticoagulation	32	32	0 (0%)	4 (12%)	2 (1, 2)
Low-dose	610	597	19 (3%)	7 (1%)	1 (1, 2)
Intermediate-dose	613	607	15 (2%)	5 (1%)	1 (1, 2)
Low-dose with aspirin	283	281	10 (4%)	5 (2%)	1 (1, 2)
Therapeutic-dose	50	50	6 (12%)	1 (2%)	1 (1, 2)
Overall	1588	1567	50 (3%)	22 (1%)	1 (1, 2)

Table 1: Summary of WHO scale at day 28 by treatment group, anticoagulation domain, FAS-ITT.

Code

p1 <- fas_itt_nona_dat %>%
  filter(!is.na(D28_who)) %>%
  group_by(CAssignment = factor(CAssignment, labels = str_replace(intervention_labels()$CAssignment, "<br>", "\n"))) %>%
  summarise(n = n()) %>%
  ggplot(., aes(CAssignment, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Anticoagulation intervention")
p2 <- fas_itt_nona_dat %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      labels = str_replace(intervention_labels()$CAssignment, "<br>", "\n")), 
    D28_who = ordered(as.integer(D28_who), levels = 1:8)
  ) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  ggplot(., aes(CAssignment, p)) +
  geom_col(aes(fill = D28_who)) +
  scale_fill_viridis_d(option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  labs(fill = "WHO scale", y = "Cumulative proportion", x = "Anticoagulation intervention")
p <- p1 | p2
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "7-2-anticoagulation.pdf"), p2, height = 3, width = 6)
p2

Figure 1: Day 28 WHO score, anticoagulation domain, FAS-ITT.

Code

p2 <- avs_itt_nona_dat %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      labels = str_replace(intervention_labels()$CAssignment, "<br>", "\n")), 
    D28_who = ordered(as.integer(D28_who), levels = 1:8)
  ) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  ggplot(., aes(CAssignment, p)) +
  geom_col(aes(fill = D28_who)) +
  scale_fill_viridis_d(option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  labs(fill = "WHO scale", y = "Cumulative proportion", x = "Anticoagulation intervention")
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "7-2-anticoagulation-avs-itt.pdf"), p2, height = 3, width = 6)
p2

Figure 2: Day 28 WHO score, anticoagulation domain, AVS-ITT.

Antiviral

There were two participants randomised to the antiviral domain who had unknown day 28 status on the WHO scale. For these two participants it was known that they were alive and out of hospital, but not whether they were grade 1 or grade 2. In the models, these participants are treated as censored.

Summary of WHO outcome by arm

save_tex_table(
  make_summary_table_antiviral(fas_itt_dat, "latex"), 
  file.path("outcomes", "secondary", "7-2-antiviral-summary"))
make_summary_table_antiviral(fas_itt_dat)

Antiviral intervention	Patients	Known	Deaths	Hospitalised	WHO, Median (Q1, Q3)
Not randomised to antiviral	1433	1414	50 (4%)	13 (1%)	1 (1, 2)
Standard of care	73	72	0 (0%)	6 (8%)	2 (1, 2)
Nafamostat	82	81	0 (0%)	3 (4%)	1 (1, 2)
Overall	1588	1567	50 (3%)	22 (1%)	1 (1, 2)

Table 2: Summary of WHO scale at day 28 by treatment group, antiviral domain, FAS-ITT.

Code

p1 <- fas_itt_nona_dat %>%
  filter(!is.na(D28_who)) %>%
  group_by(AAssignment = factor(AAssignment, labels = str_replace(intervention_labels()$AAssignment, "<br>", "\n"))) %>%
  summarise(n = n()) %>%
  ggplot(., aes(AAssignment, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Antiviral intervention")
p2 <- fas_itt_nona_dat %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    AAssignment = factor(
      AAssignment, 
      labels = str_replace(intervention_labels()$AAssignment, "<br>", "\n")), 
    D28_who = ordered(as.integer(D28_who), levels = 1:8)
  ) %>%
  group_by(AAssignment) %>%
  mutate(p = n / sum(n)) %>%
  ggplot(., aes(AAssignment, p)) +
  geom_col(aes(fill = D28_who)) +
  scale_fill_viridis_d(option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  labs(fill = "WHO scale", y = "Cumulative proportion", x = "Antiviral intervention")
p <- p1 | p2
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "7-2-antiviral.pdf"), p2, height = 3, width = 6)
p

Figure 3: Day 28 WHO score, antiviral domain, FAS-ITT.

Code

lab <- c("Usual care", "Nafamostat")
p1 <- avs_itt_nona_dat %>%
  filter(!is.na(D28_who)) %>%
  group_by(AAssignment = factor(AAssignment, labels = lab)) %>%
  summarise(n = n()) %>%
  ggplot(., aes(AAssignment, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Antiviral intervention")
p2 <- avs_itt_nona_dat %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    AAssignment = factor(AAssignment, labels = lab), 
    D28_who = ordered(as.integer(D28_who), levels = 1:8)
  ) %>%
  group_by(AAssignment) %>%
  mutate(p = n / sum(n)) %>%
  ggplot(., aes(AAssignment, p)) +
  geom_col(aes(fill = D28_who)) +
  scale_fill_viridis_d(option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  labs(fill = "WHO scale", y = "Cumulative proportion", x = "Antiviral intervention")
p <- p1 | p2
pth <- file.path("outputs", "figures", "outcomes", "secondary")
fpth <- file.path(pth, "7-2-antiviral-avs-itt.pdf") 
ggsave(fpth, p2 + coord_flip(), height = 2.5, width = 6)
system(sprintf("pdftoppm %s %s -png", fpth, gsub(".pdf", "", fpth)))
p

Figure 4: Day 28 WHO score, antiviral domain, AVS-ITT.

Code

p2 <- acs_itt_nona_dat %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    AAssignment = factor(
      AAssignment, 
      labels = str_replace(intervention_labels()$AAssignment, "<br>", "\n")), 
    D28_who = ordered(as.integer(D28_who), levels = 1:8)
  ) %>%
  group_by(AAssignment) %>%
  mutate(p = n / sum(n)) %>%
  ggplot(., aes(AAssignment, p)) +
  geom_col(aes(fill = D28_who)) +
  scale_fill_viridis_d(option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  labs(fill = "WHO scale", y = "Cumulative proportion", x = "Anticoagulation intervention")
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "7-2-antiviral-acs-itt.pdf"), p2, height = 3, width = 6)
p2

Figure 5: Day 28 WHO score, antiviral domain, ACS-ITT.

Age

Code

pdat <- all_data %>%
  filter_fas_itt() %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    agegrp = cut(AgeAtEntry, c(18, seq(25, 75, 5), 100), include.lowest = T, right = F),
    who = ordered(D28_who, levels = 1:8), 
    .drop = F) %>%
  group_by(agegrp) %>%
  mutate(p = n / sum(n))
pdat2 <- pdat %>%
  group_by(agegrp) %>%
  summarise(n = sum(n))
p1 <- ggplot(pdat2, aes(agegrp, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry") +
  geom_vline(xintercept = 8.5, linetype = 2) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45))
p2 <- ggplot(pdat, aes(agegrp, p)) +
  geom_col(aes(fill = who)) +
  labs(x = "Age", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("WHO scale", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  geom_vline(xintercept = 8.5, linetype = 2) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
  theme(legend.key.size = unit(0.5, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-2-age.pdf"), p, height = 2.5, width = 6)
p

Figure 6: Distribution of WHO scale at day 28 by age group, FAS-ITT

Sex

Code

pdat <- all_data %>%
  filter_fas_itt() %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    Sex,
    who = ordered(D28_who, levels = 1:8), 
    .drop = F) %>%
  group_by(Sex) %>%
  mutate(p = n / sum(n))
pdat2 <- pdat %>%
  group_by(Sex) %>%
  summarise(n = sum(n))
p1 <- ggplot(pdat2, aes(Sex, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  labs(y = "Number of\nparticipants")
p2 <- ggplot(pdat, aes(Sex, p)) +
  geom_col(aes(fill = who)) +
  labs(y = "Cumulative\nproportion") +
  scale_fill_viridis_d("WHO scale", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  theme(legend.key.size = unit(0.5, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-2-sex.pdf"), p, height = 2.5, width = 6)
p

Figure 7: Distribution of WHO scale at day 28 by supplemental sex requirement, FAS-ITT

Oxygen

Code

pdat <- all_data %>%
  filter_fas_itt() %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    supp_oxy = fct_explicit_na(factor(supp_oxy2, levels = 0:1, labels = c("Not required", "Required"))),
    who = ordered(D28_who, levels = 1:8), 
    .drop = F) %>%
  group_by(supp_oxy) %>%
  mutate(p = n / sum(n))
pdat2 <- pdat %>%
  group_by(supp_oxy) %>%
  summarise(n = sum(n))
p1 <- ggplot(pdat2, aes(supp_oxy, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  labs(y = "Number of\nparticipants",
       x = "Supplemental oxygen")
p2 <- ggplot(pdat, aes(supp_oxy, p)) +
  geom_col(aes(fill = who)) +
  labs(x = "Supplemental oxygen", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("WHO scale", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  theme(legend.key.size = unit(0.5, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-2-oxygen.pdf"), p, height = 2.5, width = 6)
p

Figure 8: Distribution of WHO scale at day 28 by supplemental oxygen requirement, FAS-ITT

Country

Code

pdat <- all_data %>%
  filter_fas_itt() %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand"),
                     labels = c("India", "Australia", "Nepal", "New\nZealand")),
    who = ordered(D28_who, levels = 1:8), 
    .drop = F) %>%
  group_by(Country) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
pdat2 <- pdat %>%
  group_by(Country) %>%
  summarise(n = sum(n))
p1 <- ggplot(pdat2, aes(Country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")
p2 <- ggplot(pdat, aes(Country, p)) +
  geom_col(aes(fill = who)) +
  labs(x = "Country", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("WHO scale", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  theme(legend.key.size = unit(0.5, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-2-country.pdf"), p, height = 2.5, width = 6)
p

Figure 9: Distribution of WHO scale at day 28 by country, FAS-ITT.

Site

Code

pdat <- all_data %>%
  filter_fas_itt() %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand"),
                     labels = c("India", "Australia", "Nepal", "New\nZealand")),
    Site = fct_infreq(Location),
    who = ordered(D28_who, levels = 1:8)) %>%
  complete(who = ordered(1:8), nesting(Country, Site), fill = list(n = 0)) %>%
  group_by(Country, Site) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup() %>%
  mutate(
    Country = droplevels(Country),
    Site = droplevels(Site)
  )
pdat2 <- pdat %>%
  group_by(Country, Site) %>%
  summarise(n = sum(n)) %>%
  ungroup()
p1 <- ggplot(pdat2, aes(Site, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(pdat, aes(Site, p)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col(aes(fill = who)) +
  labs(x = "Study site", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("WHO scale", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F, ncol = 1)) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
  theme(legend.key.size = unit(0.25, "lines"))
p <- p1 / p2 +
  plot_layout(guides = 'collect')
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-2-country-site.pdf"), p, height = 4, width = 6.25)
p

Figure 10: Distribution of WHO scale by study site, FAS-ITT.

Calendar Time

Code

pdat <- fas_itt_nona_dat %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    yr = year(RandDate), mth = month(RandDate),
    who = ordered(D28_who, levels = 1:8), 
    .drop = F) %>%
  group_by(yr, mth) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p1 <- pdat %>%
  group_by(yr, mth) %>%
  summarise(n = sum(n)) %>%
  ggplot(., aes(mth, n))  +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p2 <- ggplot(pdat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
  geom_col(aes(fill = who)) +
  labs(x = "Calendar date (month of year)", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("WHO scale", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  theme(legend.key.size = unit(0.5, "lines")) +
  scale_x_continuous(breaks = 1:12)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-2-calendar-time.pdf"), p, height = 2, width = 6)
p

Figure 11: Distribution of WHO scale by calendar time, FAS-ITT.

Sample Cumulative Logits

Proportional odds looks reasonable for all logits except 1.

Plot sample cumulative logits

trt_counts <- fas_itt_nona_dat %>%
  dplyr::count(CAssignment, D28_who) %>%
  complete(CAssignment, D28_who, fill = list(n = 0)) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n))
trt_logit <- trt_counts %>% 
  group_by(CAssignment) %>% 
  mutate(clogit = logit(cumsum(p))) %>%
  group_by(D28_who) %>%
  mutate(rel_clogit = clogit - mean(clogit)) %>%
  filter(D28_who != 28)
ggplot(trt_logit, aes(D28_who, rel_clogit)) +
  facet_wrap( ~ CAssignment) +
  geom_point() +
  labs(y = "Relative (to mean) sample cumulative logit")

Inspect sample cumulative logits, anticoagulation.

Plot sample cumulative logits

trt_counts <- fas_itt_nona_dat %>%
  dplyr::count(AAssignment, D28_who) %>%
  complete(AAssignment, D28_who, fill = list(n = 0)) %>%
  group_by(AAssignment) %>%
  mutate(p = n / sum(n))
trt_logit <- trt_counts %>% 
  group_by(AAssignment) %>% 
  mutate(clogit = logit(cumsum(p))) %>%
  group_by(D28_who) %>%
  mutate(rel_clogit = clogit - mean(clogit)) %>%
  filter(D28_who != 28)
ggplot(trt_logit, aes(D28_who, rel_clogit)) +
  facet_wrap( ~ AAssignment) +
  geom_point() +
  labs(y = "Relative (to mean) sample cumulative logit")

Inspect sample cumulative logits, antiviral.

Modelling

FAS-ITT

Fit primary model

res <- fit_primary_model(
  fas_itt_nona_dat,
  ordmod,
  outcome = "D28_who",
  intercept = FALSE
)
# res2 <- fit_primary_model2(
#   fas_itt_dat |> filter(!is.na(D28_who_low)),
#   ordmodcens,
#   outcome = c("D28_who_low", "D28_who_high"),
#   intercept = FALSE
# )
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Ineligible aspirin", "Age \u2265 60", "Female", "Oxygen requirement", "Australia/New Zealand", "Nepal")

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-2-primary-model-fas-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Nafamostat	0.67	(0.35, 1.29)	0.71 (0.24)	0.88
Intermediate-dose	0.80	(0.60, 1.07)	0.81 (0.12)	0.94
Low-dose with aspirin	0.73	(0.51, 1.05)	0.75 (0.14)	0.95
Therapeutic-dose	1.70	(0.84, 3.41)	1.81 (0.67)	0.07
Ineligible aspirin	1.64	(0.73, 3.51)	1.76 (0.72)	0.11
Age ≥ 60	2.50	(1.91, 3.27)	2.52 (0.35)	0.00
Female	0.94	(0.73, 1.22)	0.95 (0.12)	0.67
Oxygen requirement	2.11	(1.59, 2.82)	2.14 (0.32)	0.00
Australia/New Zealand	1.30	(0.42, 3.91)	1.52 (0.93)	0.32
Nepal	0.55	(0.17, 1.96)	0.68 (0.51)	0.84

Posterior summaries for model parameters (fixed-effects), FAS-ITT.

Odds ratio posterior densities

p <- plot_or_densities(c(res$drws$AOR, res$drws$COR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "primary", "7-2-primary-model-fas-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Figure 12: Posterior densities for odds ratio contrasts.

Odds ratio posterior summary for epoch and site

p <- plot_epoch_site_terms(
  res$drws$gamma_epoch,
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-2-primary-model-epoch-site-terms-fas-itt.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Summary of epoch and site posterior odds ratios.

Posterior Predictive

Posterior predictive figure

y_ppc <- res$drws$y_ppc
ppc_dat <- bind_cols(fas_itt_nona_dat, tibble(y_ppc = y_ppc))
grp_ppc1 <- function(grp) {
  ppc_dat  %>%
  group_by(grp = !!grp) %>%
  summarise(
    y_1 = mean(D28_who == 1),
    ypp_1 = rvar_mean(y_ppc == 1),
    y_2 = mean(D28_who <= 2),
    ypp_2 = rvar_mean(y_ppc <= 2),
    y_3 = mean(D28_who <= 3),
    ypp_3 = rvar_mean(y_ppc <= 3),
    y_4 = mean(D28_who <= 4),
    ypp_4 = rvar_mean(y_ppc <= 4),
    y_5 = mean(D28_who <= 5),
    ypp_5 = rvar_mean(y_ppc <= 5),
    y_6 = mean(D28_who <= 6),
    ypp_6 = rvar_mean(y_ppc <= 6),
    y_7 = mean(D28_who <= 7),
    ypp_7 = rvar_mean(y_ppc <= 7),
    y_8 = mean(D28_who <= 8),
    ypp_8 = rvar_mean(y_ppc <= 8)
  ) %>%
  pivot_longer(y_1:ypp_8, names_to = c("response", "who"), names_sep = "_", values_to = "posterior") %>%
  mutate(who = as.numeric(who))
}
grp_ppc2 <- function(grp) {
  ppc_dat  %>%
  group_by(grp = !!grp) %>%
  summarise(
    y_1 = mean(D28_who == 1),
    ypp_1 = rvar_mean(y_ppc == 1),
    y_2 = mean(D28_who == 2),
    ypp_2 = rvar_mean(y_ppc == 2),
    y_3 = mean(D28_who == 3),
    ypp_3 = rvar_mean(y_ppc == 3),
    y_4 = mean(D28_who == 4),
    ypp_4 = rvar_mean(y_ppc == 4),
    y_5 = mean(D28_who == 5),
    ypp_5 = rvar_mean(y_ppc == 5),
    y_6 = mean(D28_who == 6),
    ypp_6 = rvar_mean(y_ppc == 6),
    y_7 = mean(D28_who == 7),
    ypp_7 = rvar_mean(y_ppc == 7),
    y_8 = mean(D28_who == 8),
    ypp_8 = rvar_mean(y_ppc == 8)
  ) %>%
  pivot_longer(y_1:ypp_8, names_to = c("response", "who"), 
               names_sep = "_", values_to = "posterior") %>%
  mutate(who = as.numeric(who))
}
plot_grp_ppc <- function(dat) {
  ggplot(dat %>% filter(response == "ypp"), aes(x = who)) +
    facet_wrap( ~ grp, nrow = 1) +
    stat_slabinterval(aes(ydist = posterior))  +
    geom_point(data = dat %>% filter(response == "y"), 
               aes(x = who, y = mean(posterior)),
               colour = "red",
               shape = 23) +
    labs(x = "WHO outcome", 
         y = "Probability")
}
plot_grp_ppc <- function(dat, lab = "", xlab = "Probability") {
  ggplot(dat %>% filter(response == "ypp"), aes(y = who)) +
    facet_wrap( ~ grp, nrow = 1) +
    stat_interval(aes(xdist = posterior), size = 2)  +
    geom_point(data = dat %>% filter(response == "y"), 
               aes(y = who, x = mean(posterior)),
               colour = "red",
               shape = 23) +
    scale_x_continuous(xlab, breaks = c(0, 0.5),
                       sec.axis = sec_axis(~ . , name = lab, breaks = NULL, labels = NULL)) +
    scale_y_continuous("WHO\noutcome", breaks = c(1,3,5,7)) +
    theme(strip.text = element_text(size = rel(0.7)),
          axis.title.x = element_text(size = rel(0.7)),
          axis.text.x = element_text(size = rel(0.65)),
          axis.title.y = element_text(size = rel(0.75)),
          axis.title.x.bottom = element_blank())
}
pp_epoch <- grp_ppc2(sym("epoch")) %>% 
  mutate(grp = fct_inorder(factor(grp)))
pp_A <- grp_ppc2(sym("AAssignment"))
pp_C <- grp_ppc2(sym("CAssignment"))
pp_ctry <- grp_ppc2(sym("Country"))
pp_site <- grp_ppc2(sym("site")) %>%
  left_join(ppc_dat %>% dplyr::count(site, Country), by = c("grp" = "site"))
p0 <- plot_grp_ppc(pp_A, "Antiviral", "")
p1 <- plot_grp_ppc(pp_C, "Anticoagulation", "")
p2 <- plot_grp_ppc(pp_ctry, "Country", "") 
p3 <- plot_grp_ppc(pp_epoch, "Epoch", "")
p4 <- plot_grp_ppc(pp_site %>% filter(Country == "IN"), "Sites India", "")
p5 <- plot_grp_ppc(pp_site %>% filter(Country == "AU"), "Sites Australia", "")
p6 <- plot_grp_ppc(pp_site %>% filter(Country == "NP"), "Sites Nepal", "")
p7 <- plot_grp_ppc(pp_site %>% filter(Country == "NZ"), "Sites New Zealand", "")
p <- (p0 | p1 | p2) / p3 / p4 / p5 / (p6 | p7)+
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")
pth <- file.path("outputs", "figures", "outcomes", "secondary",
                 "7-2-primary-model-fas-itt-ppc.pdf")
ggsave(pth, p, width = 6, height = 5.75)
p

Figure 13: Posterior predictive check, FAS-ITT.

ACS-ITT

Fit primary model

res <- fit_primary_model(
  acs_itt_nona_dat,
  ordmod,
  outcome = "D28_who",
  intercept = FALSE
)
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Ineligible aspirin", "Age \u2265 60", "Female", "Oxygen requirement", "Australia/New Zealand", "Nepal")

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-2-primary-model-acs-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Nafamostat	0.66	(0.31, 1.41)	0.71 (0.28)	0.86
Intermediate-dose	0.80	(0.61, 1.07)	0.81 (0.12)	0.93
Low-dose with aspirin	0.73	(0.51, 1.06)	0.75 (0.14)	0.95
Therapeutic-dose	1.72	(0.84, 3.43)	1.82 (0.67)	0.07
Ineligible aspirin	1.67	(0.76, 3.63)	1.81 (0.75)	0.10
Age ≥ 60	2.36	(1.80, 3.10)	2.38 (0.33)	0.00
Female	0.94	(0.73, 1.21)	0.95 (0.12)	0.68
Oxygen requirement	2.21	(1.65, 2.95)	2.23 (0.33)	0.00
Australia/New Zealand	1.14	(0.34, 3.54)	1.35 (0.86)	0.41
Nepal	0.59	(0.18, 2.09)	0.73 (0.54)	0.81

Posterior summaries for model parameters (fixed-effects), ACS-ITT.

Odds ratio posterior densities

p <- plot_or_densities(c(res$drws$AOR, res$drws$COR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "primary", "7-2-primary-model-acs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Figure 14: Posterior densities for odds ratio contrasts, ACS-ITT.

Odds ratio posterior summary for epoch and site

p <- plot_epoch_site_terms(
  res$drws$gamma_epoch,
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-2-primary-model-epoch-site-terms-acs-itt.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Figure 15: Summary of epoch and site posterior odds ratios, ACS-ITT.

AVS-ITT

Code

p1 <- avs_itt_nona_dat %>%
  filter(!is.na(D28_who)) %>%
  group_by(CAssignment = factor(CAssignment, labels = str_replace(intervention_labels()$CAssignment, "<br>", "\n"))) %>%
  summarise(n = n()) %>%
  ggplot(., aes(CAssignment, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Anticoagulation intervention")
p2 <- fas_itt_nona_dat %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      labels = str_replace(intervention_labels()$CAssignment, "<br>", "\n")), 
    D28_who = ordered(as.integer(D28_who), levels = 1:8)
  ) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  ggplot(., aes(CAssignment, p)) +
  geom_col(aes(fill = D28_who)) +
  scale_fill_viridis_d(option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  labs(fill = "WHO scale", y = "Cumulative proportion", x = "Anticoagulation intervention")
p <- p1 | p2
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "7-2-anticoagulation-avs-itt.pdf"), p2, height = 3, width = 6)
p2

Pre-specified model

• excludes epoch

Fit primary model - AVS-ITT

res <- fit_primary_model(
  avs_itt_nona_dat,
  ordmod_site,
  outcome = "D28_who",
  vars = c("agegte60", "sexF", "supp_oxy2", "crp_tertile", "ctry"),
  beta_sd_var = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 1),
  intercept = FALSE
)
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Age \u2265 60", "Female", "Required oxygen", "CRP (2nd tertile)", "CRP (3rd tertile)", "CRP (unknown)", "Nepal")

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table_rev(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-2-primary-model-avs-itt-summary-table-pre-spec")
odds_ratio_summary_table_rev(c(res$drws$AOR, res$drws$COR, res$drws$OR))

Parameter	Median	95% CrI	Mean (SD)	Pr(OR > 1)
Nafamostat	0.64	(0.30, 1.35)	0.69 (0.28)	0.12
Intermediate-dose	1.21	(0.52, 2.82)	1.33 (0.61)	0.67
Low-dose with aspirin	1.13	(0.25, 4.77)	1.48 (1.26)	0.56
Therapeutic-dose	1.97	(0.61, 6.42)	2.36 (1.57)	0.87
Age ≥ 60	4.06	(1.82, 9.35)	4.45 (1.96)	1.00
Female	1.27	(0.61, 2.59)	1.36 (0.51)	0.74
Required oxygen	0.95	(0.42, 2.16)	1.03 (0.45)	0.45
CRP (2nd tertile)	0.85	(0.33, 2.09)	0.94 (0.46)	0.36
CRP (3rd tertile)	0.95	(0.38, 2.35)	1.06 (0.52)	0.46
CRP (unknown)	0.96	(0.26, 3.48)	1.19 (0.86)	0.48
Nepal	0.40	(0.08, 2.21)	0.59 (0.65)	0.14

Posterior summaries for model parameters (fixed-effects), AVS-ITT.

Code

p <- plot_or_densities(c(res$drws$AOR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "primary", "7-2-primary-model-avs-itt-odds-ratio-densities-pre-spec.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Figure 16: Posterior densities for odds ratio contrasts, AVS-ITT.

Odds ratio summary for epoch and site

p <- plot_site_terms(
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-2-primary-model-site-terms-avs-itt-pre-spec.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Summary of epoch and site posterior odds ratios.

Diagnostics

Code

res$fit$summary(c("beta", "gamma_site", "tau_site")) |> print(n = Inf)

# A tibble: 28 × 10
   variable     mean   median    sd   mad      q5    q95  rhat ess_bulk ess_tail
   <chr>       <num>    <num> <num> <num>   <num>  <num> <num>    <num>    <num>
 1 beta[1]  -0.313   -0.312   0.273 0.273 -0.766   0.133  1.00   26668.   15387.
 2 beta[2]   0.270    0.269   0.484 0.484 -0.530   1.07   1.00   20445.   16197.
 3 beta[3]   0.404    0.400   0.606 0.609 -0.591   1.40   1.00   24642.   15748.
 4 beta[4]  -0.283   -0.282   0.375 0.372 -0.904   0.336  1.00   22689.   16019.
 5 beta[5]  -0.151   -0.156   0.438 0.435 -0.868   0.577  1.00   20373.   16485.
 6 beta[6]   1.41     1.40    0.417 0.416  0.724   2.09   1.00   26578.   16205.
 7 beta[7]   0.238    0.240   0.368 0.368 -0.371   0.840  1.00   28203.   14959.
 8 beta[8]  -0.0535  -0.0526  0.416 0.416 -0.736   0.633  1.00   21077.   15124.
 9 beta[9]  -0.168   -0.168   0.463 0.461 -0.935   0.600  1.00   18615.   16712.
10 beta[10] -0.0463  -0.0463  0.466 0.462 -0.826   0.714  1.00   18432.   15806.
11 beta[11] -0.0429  -0.0387  0.662 0.662 -1.15    1.04   1.00   18926.   15337.
12 beta[12] -0.903   -0.909   0.843 0.827 -2.26    0.506  1.00   22678.   14762.
13 gamma_s…  0.290    0.274   0.531 0.507 -0.557   1.18   1.00   23231.   16993.
14 gamma_s… -0.552   -0.514   0.579 0.570 -1.56    0.324  1.00   19077.   15993.
15 gamma_s…  0.270    0.253   0.568 0.536 -0.642   1.23   1.00   25020.   16866.
16 gamma_s…  0.683    0.651   0.566 0.557 -0.176   1.66   1.00   18662.   16827.
17 gamma_s…  0.00744  0.00935 0.660 0.607 -1.08    1.09   1.00   27440.   16332.
18 gamma_s… -0.197   -0.182   0.495 0.477 -1.03    0.584  1.00   21609.   15179.
19 gamma_s…  0.580    0.544   0.609 0.599 -0.348   1.63   1.00   23225.   16032.
20 gamma_s…  0.700    0.660   0.627 0.609 -0.241   1.79   1.00   19815.   16298.
21 gamma_s… -0.365   -0.323   0.632 0.598 -1.46    0.594  1.00   23882.   16202.
22 gamma_s… -0.100   -0.0898  0.572 0.540 -1.06    0.828  1.00   24378.   17150.
23 gamma_s… -1.15    -1.13    0.491 0.486 -2.00   -0.379  1.00   13626.    9472.
24 gamma_s… -0.483   -0.437   0.593 0.572 -1.53    0.414  1.00   19914.   15624.
25 gamma_s… -0.393   -0.146   0.944 0.565 -2.24    0.789  1.00   17506.   17383.
26 gamma_s… -0.636   -0.214   1.40  0.612 -3.21    0.691  1.00   18226.   16951.
27 tau_sit…  0.814    0.783   0.298 0.270  0.393   1.35   1.00    7976.    7923.
28 tau_sit…  1.00     0.724   1.01  0.682  0.0609  2.84   1.00   12766.    8421.

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8817 0.8833 0.8452 0.9155 0.8299 0.9196 0.8230 0.8229

Trace plots

Code

mcmc_trace(res$drws["beta"])

Code

mcmc_trace(res$drws["gamma_site"])

Code

mcmc_trace(res$drws["tau_site"])

Reduced Model

• excludes epoch, region, and site

Code

res <- ASCOTr::fit_primary_model(
  avs_itt_nona_dat,
  ordmod0,
  outcome = "D28_who",
  vars = c("agegte60", "sexF", "supp_oxy2", "crp_tertile"),
  beta_sd_var = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5),
  intercept = FALSE
)
# res2 <- fit_primary_model2(
#   avs_itt_dat |> filter(!is.na(D28_who_low)),
#   ordmodcens0,
#   outcome = c("D28_who_low", "D28_who_high"),
#   vars = c("agegte60", "sexF", "supp_oxy2", "crp_tertile"),
#   beta_sd_var = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5),
#   intercept = FALSE
# )
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Age \u2265 60", "Female", "Oxygen requirement", "CRP (2nd tertile)", "CRP (3rd tertile)", "CRP (unknown)")

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-2-primary-model-avs-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Nafamostat	0.72	(0.36, 1.45)	0.77 (0.28)	0.82
Intermediate-dose	1.27	(0.56, 2.87)	1.38 (0.60)	0.28
Low-dose with aspirin	0.80	(0.19, 3.20)	1.03 (0.83)	0.62
Therapeutic-dose	1.97	(0.65, 5.87)	2.30 (1.40)	0.11
Age ≥ 60	3.77	(1.77, 8.24)	4.08 (1.70)	0.00
Female	1.29	(0.65, 2.53)	1.37 (0.49)	0.23
Oxygen requirement	1.16	(0.55, 2.54)	1.27 (0.51)	0.34
CRP (2nd tertile)	0.81	(0.35, 1.86)	0.89 (0.39)	0.69
CRP (3rd tertile)	0.79	(0.33, 1.84)	0.86 (0.39)	0.71
CRP (unknown)	0.84	(0.26, 2.54)	0.98 (0.61)	0.62

Posterior summaries for model parameters (fixed-effects), AVS-ITT.

Odds ratio posterior densities

p <- plot_or_densities(c(res$drws$AOR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "primary", "7-2-primary-model-avs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Figure 17: Posterior densities for odds ratio contrasts, AVS-ITT.

Diagnostics

Code

res$fit$summary(c("beta"))

# A tibble: 11 × 10
   variable    mean  median    sd   mad     q5   q95  rhat ess_bulk ess_tail
   <chr>      <num>   <num> <num> <num>  <num> <num> <num>    <num>    <num>
 1 beta[1]  -0.233  -0.232  0.251 0.252 -0.643 0.178  1.00   19316.   15028.
 2 beta[2]   0.427   0.425  0.411 0.411 -0.246 1.10   1.00   18249.   15413.
 3 beta[3]   0.641   0.635  0.571 0.567 -0.287 1.59   1.00   20616.   15338.
 4 beta[4]  -0.199  -0.195  0.358 0.356 -0.784 0.394  1.00   18357.   15409.
 5 beta[5]   0.0320  0.0306 0.412 0.409 -0.642 0.708  1.00   17076.   15033.
 6 beta[6]   1.33    1.33   0.394 0.392  0.687 1.98   1.00   21643.   15483.
 7 beta[7]   0.255   0.258  0.346 0.345 -0.322 0.826  1.00   21612.   14448.
 8 beta[8]   0.160   0.151  0.388 0.389 -0.473 0.808  1.00   17776.   14276.
 9 beta[9]  -0.210  -0.209  0.426 0.426 -0.919 0.490  1.00   16159.   15271.
10 beta[10] -0.244  -0.242  0.432 0.430 -0.954 0.470  1.00   16257.   15053.
11 beta[11] -0.186  -0.173  0.583 0.584 -1.16  0.755  1.00   19310.   15520.

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0188 0.8839 0.9600 0.9121 0.9629 0.9493 0.9682 0.9562

Trace plots

Code

mcmc_trace(res$drws["beta"])

Country

Rather than group Australia and New Zealand, keep them as separate countries in the model.

• excludes epoch

Fit primary model - AVS-ITT

res <- fit_primary_model(
  avs_itt_nona_dat,
  ordmod_site,
  outcome = "D28_who",
  vars = c("agegte60", "sexF", "supp_oxy2", "crp_tertile", "ctry2"),
  beta_sd_var = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 1, 1),
  intercept = FALSE,
  ctry_var = "ctry2_num",
  site_var = "site2_num"
)
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Age \u2265 60", "Female", "Required oxygen", "CRP (2nd tertile)", "CRP (3rd tertile)", "CRP (unknown)", "Nepal", "New Zealand")

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table_rev(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-2-primary-model-avs-itt-summary-table-ctry2")
odds_ratio_summary_table_rev(c(res$drws$AOR, res$drws$COR, res$drws$OR))

Parameter	Median	95% CrI	Mean (SD)	Pr(OR > 1)
Nafamostat	0.65	(0.30, 1.37)	0.70 (0.28)	0.13
Intermediate-dose	1.20	(0.50, 2.91)	1.32 (0.63)	0.66
Low-dose with aspirin	1.11	(0.24, 4.90)	1.47 (1.27)	0.55
Therapeutic-dose	1.98	(0.60, 6.61)	2.39 (1.67)	0.87
Age ≥ 60	4.01	(1.80, 9.34)	4.40 (1.97)	1.00
Female	1.25	(0.60, 2.60)	1.35 (0.52)	0.73
Required oxygen	0.94	(0.41, 2.15)	1.03 (0.46)	0.44
CRP (2nd tertile)	0.84	(0.34, 2.10)	0.94 (0.47)	0.35
CRP (3rd tertile)	0.96	(0.38, 2.42)	1.07 (0.53)	0.46
CRP (unknown)	0.92	(0.24, 3.34)	1.13 (0.83)	0.45
Nepal	0.40	(0.08, 2.19)	0.59 (0.63)	0.14
New Zealand	0.70	(0.18, 2.76)	0.89 (0.72)	0.29

Posterior summaries for model parameters (fixed-effects), AVS-ITT.

Code

p <- plot_or_densities(c(res$drws$AOR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "primary", "7-2-primary-model-avs-itt-odds-ratio-densities-ctry2.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Figure 18: Posterior densities for odds ratio contrasts, AVS-ITT.

Odds ratio summary for epoch and site

p <- plot_site_terms(
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("Australia", "Nepal", "New Zealand"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-2-primary-model-site-terms-avs-itt-ctry2.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Summary of epoch and site posterior odds ratios.

Diagnostics

Code

res$fit$summary(c("beta", "gamma_site", "tau_site")) |> print(n = Inf)

# A tibble: 31 × 10
   variable     mean   median    sd   mad      q5    q95  rhat ess_bulk ess_tail
   <chr>       <num>    <num> <num> <num>   <num>  <num> <num>    <num>    <num>
 1 beta[1]  -0.310   -0.308   0.275 0.274 -0.768   0.138  1.00   23858.   15565.
 2 beta[2]   0.267    0.269   0.489 0.492 -0.541   1.06   1.00   18212.   14603.
 3 beta[3]   0.419    0.417   0.616 0.615 -0.589   1.43   1.00   21616.   15339.
 4 beta[4]  -0.276   -0.275   0.384 0.384 -0.912   0.358  1.00   20330.   16110.
 5 beta[5]  -0.156   -0.156   0.444 0.442 -0.882   0.575  1.00   18001.   15759.
 6 beta[6]   1.39     1.39    0.421 0.422  0.703   2.09   1.00   25062.   14679.
 7 beta[7]   0.227    0.226   0.374 0.370 -0.391   0.842  1.00   24377.   15882.
 8 beta[8]  -0.0598  -0.0647  0.423 0.425 -0.750   0.644  1.00   19573.   15236.
 9 beta[9]  -0.173   -0.169   0.463 0.462 -0.940   0.596  1.00   17597.   15705.
10 beta[10] -0.0436  -0.0451  0.468 0.460 -0.806   0.729  1.00   18346.   16381.
11 beta[11] -0.0929  -0.0887  0.662 0.660 -1.19    0.987  1.00   17953.   15401.
12 beta[12] -0.895   -0.906   0.837 0.829 -2.25    0.492  1.00   22065.   16660.
13 beta[13] -0.359   -0.359   0.688 0.666 -1.47    0.770  1.00   16744.   14598.
14 gamma_s…  0.294    0.280   0.541 0.525 -0.565   1.21   1.00   21222.   16707.
15 gamma_s… -0.453   -0.422   0.599 0.579 -1.48    0.469  1.00   19625.   15392.
16 gamma_s…  0.276    0.263   0.585 0.563 -0.652   1.26   1.00   22732.   16235.
17 gamma_s…  0.715    0.683   0.581 0.573 -0.184   1.72   1.00   17866.   17037.
18 gamma_s…  0.00234  0.0114  0.679 0.628 -1.12    1.09   1.00   24160.   16584.
19 gamma_s…  0.611    0.580   0.623 0.603 -0.354   1.68   1.00   22447.   16294.
20 gamma_s…  0.738    0.700   0.643 0.626 -0.234   1.84   1.00   20159.   16639.
21 gamma_s… -0.404   -0.363   0.673 0.640 -1.57    0.617  1.00   18941.   15848.
22 gamma_s… -0.116   -0.108   0.590 0.556 -1.09    0.840  1.00   21138.   16256.
23 gamma_s… -1.19    -1.17    0.492 0.487 -2.04   -0.424  1.00   13748.   12864.
24 gamma_s… -0.514   -0.478   0.602 0.580 -1.54    0.413  1.00   18226.   15940.
25 gamma_s… -0.400   -0.158   0.943 0.559 -2.24    0.795  1.00   17198.   17892.
26 gamma_s… -0.652   -0.217   1.53  0.613 -3.24    0.708  1.00   17438.   17560.
27 gamma_s…  0.0298   0.00912 0.620 0.410 -0.986   1.09   1.00   16004.   15925.
28 gamma_s… -0.424   -0.125   1.14  0.498 -2.46    0.752  1.00   19652.   16603.
29 tau_sit…  0.863    0.825   0.313 0.286  0.429   1.42   1.00    9596.   10382.
30 tau_sit…  1.02     0.721   1.14  0.684  0.0605  2.91   1.00   14039.    9197.
31 tau_sit…  0.819    0.604   0.806 0.553  0.0568  2.28   1.00   12698.    9694.

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8756 0.9545 0.8733 0.8431 0.9059 0.8833 0.8768 0.9305

Trace plots

Code

mcmc_trace(res$drws["beta"])

Code

mcmc_trace(res$drws["gamma_site"])

Code

mcmc_trace(res$drws["tau_site"])